require(mvtnorm)
require(OpenMx)
require(psych)
require(MASS)

rAE <- 0                                                 # Set the level of the correlation betwween the A and C Varaince Components
nsim <- 200000                                           # Set the number of simulations (We used 200000)

MZdatasets <- list(NULL)
DZdatasets <- list(NULL)

ests <- matrix(NA, nsim, 6)

nv <- 1                                                     # Set the number of variables that you will be using (per person)
Vars <- c('t')                                              # Name the variable
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")        # Create an object that will select the variables
ntv <- nv*2                                                 # Set the total number of variable that you will use (1 var * 2 twins)




mzmat <- matrix(c(  1,  0,rAE,rAE,
                    0,  1,  0,  0,
                  rAE,  0,  1,  0,
                  rAE,  0,  0,  1),4)                       # Create a matrix of correlations for the genetic and environmental variance components for each twin (MZ)
				                                            #  Note that because the "Genes" the same there is only 1 factor score estimated for each twin pair
				                                            #  Note that because the "Shared Environment" the same there is only 1 factor score estimated for each twin pair
															
dzmat <- matrix(c(  1, .5,  0,rAE,  0,  
                   .5,  1,  0,  0,rAE,
                    0,  0,  1,  0,  0,
                  rAE,  0,  0,  1,  0,
                    0,rAE,  0,  0,  1),5)                   # Create a matrix of correlations for the genetic and environmental variance components for each twin (MZ)
				                                            #  Note that because the "Shared Environment" the same there is only 1 factor score estimated for each twin pair



 for(sim in 1:nsim) {

	 mzSim  <- mvrnorm(n=1000, c(0,0,0,0), mzmat , empirical=T)        # Simulate genetic and environmental factor scores for each MZ twin
	 dzSim  <- mvrnorm(n=1000, c(0,0,0,0,0), dzmat , empirical=T)      # Simulate genetic and environmental factor scores for each DZ twin

	 VarCompE <- runif(1,.2,.5)                                            # Set the Genetic Variance Component
	 VarCompA <- runif(1,0,(1-VarCompE))                                   # Set the Common Environmental Variance Component
	 VarCompC <- 1-VarCompA -VarCompE                                      # Set the Unique Environmental Variance Component


	 aPath <- sqrt(VarCompA)
	 cPath <- sqrt(VarCompC)
	 ePath <- sqrt(VarCompE)

	 MZdata <- cbind((mzSim[,1]*aPath+
	                  mzSim[,2]*cPath+
	                  mzSim[,3]*ePath),
	                 (mzSim[,1]*aPath+
	                  mzSim[,2]*cPath+
	                  mzSim[,4]*ePath))

	 DZdata <- cbind((dzSim[,1]*aPath+
	                  dzSim[,3]*cPath+
	                  dzSim[,4]*ePath),
	                 (dzSim[,2]*aPath+
	                  dzSim[,3]*cPath+
	                  dzSim[,5]*ePath))

	 colnames (MZdata) <- selVars
	 colnames (DZdata) <- selVars


	 colnames (MZdata) <- selVars
	 colnames (DZdata) <- selVars
		
		MZdatasets[[sim]] <- MZdata	                                           # Save the MZ data
		DZdatasets[[sim]] <- DZdata	                                           # Save the DZ data	
		ests[sim,1:6] <- c(aPath,cPath,ePath,VarCompA,VarCompC,VarCompE)       # Save the simulated parameter values for each rep
		
	}
colnames(ests) <- c("aPath","cPath","ePath","VarCompA","VarCompC","VarCompE")   

#MZdatasets
#DZdatasets